function K = commute(n, m)
% Finds the commutation matrix K(m,n) as in Magnus and Neudecker (2002)
% such that
%    K(m,n)*vec(A) = vec(A')

K = reshape(kron(vec(eye(n)), eye(m)), n*m, n*m);
end
